Interface dynamics from experimental data 
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A novel algorithm is envisaged to extract the coupling 
parameters of the Kardar-Parisi-Zhang (KPZ) equation from 
experimental data. The method hinges on the Fokker-Planck 
equation combined with a classical least-square error proce- 
dure. It takes properly into account the fluctuations of surface 
height through a deterministic equation for space correlations. 
We apply it to the 1+1 KPZ equation and carefully compare 
its results with those obtained by previous investigations. Un- 
like previous approaches, our method does not require large 
sizes and is stable under a modification of sampling time of 
observations. Shortcomings associated to standard discretiza- 
tions of the continuous KPZ equation are also pointed out and 
possible future perspectives are finally analyzed. 

64.60.Ht,05.40.+j,05.70.Ln 



I. INTRODUCTION 

Inverse techniques have a wide range of applicabil- 
ity ranging from geophysics to nonlinear time analysis 
and statistics jlj. The common philosophy behind these 
methods is the extraction of equations of motion start- 
ing from successive experimental time series of some dy- 
namical variable in addition to basic assumptions such as 
determinism. If a reasonably general form of the equa- 
tions is guessed either by symmetry arguments or by gen- 
eral considerations, the "true" parameters are then de- 
termined by minimizing a cost function quantifying the 
distance between experimental observations and corre- 
sponding reconstructed quantities, the latter being im- 
plicitly dependent upon the parameters. Among such 
approaches, the Least Square method is the most popu- 
lar one. 

A typical system which can be treated using recon- 
struction techniques is the case where an observational 
noise is superimposed to a standard deterministic evolu- 
tion. In this case the system is expected to evolve un- 
der the action of deterministic system and stochasticity 
come only from our measurement apparatus. The par- 
ticular case where the dynamics underlying the system 
is chaotic has also received considerable attention due to 
its widespread occurence in natural systems Q , and the 
importance of treating the presence of the noise with due 
care has already been emphasized (3) . 

The alternative possibility of dynamical noise occurs 
whenever the noise is a built-in component of the equa- 
tions of motion. This is a far more difficult problem since 



one has to deal with stochastic rather than deterministic 
equations. An important such case, which is rather ubiq- 
uitous in nature, is the Langevin dynamics where vari- 
ables evolve subject both to dissipative generalized forces 
and to a fluctuating part Q. In this latter instance, the 
presence of dynamical noise can drastically modify the 
dynamics and hence hampers the efficiency of the usual 
reconstruction techniques based on deterministic ideas 

1- 

In our work, we focus on a particular class of Langevin 
dynamics which has its origin in a seminal paper on inter- 
face dynamics |(| but has ever since displayed relations to 
a variety of physical systems such as for instance bacte- 
rial colonial growth, immiscible fluids, directed polymer 
and superconductors 0. 

The Kardar-Parisi-Zhang (KPZ) equation || was in- 
troduced as a coarse-grained mesoscopic description for 
the growth of a rough surface under the deposition of 
particles driven by gravity. The crucial ingredient intro- 
duced in the KPZ and not present in the corresponding 
linear counterpart, namely the Edward- Wilkinson (EW) 
equation j^] , is a non- linear term which takes into account 
the fact that the growth is normal to the surface. The 
KPZ equation can be mapped into various other models. 
A Cole-Hopf change of variables maps it into a directed 
polymer diffusion equation subject to a random potential 
while the identification of the local gradient with a 
velocity leads to the Burgers's equation for a vorticity- 
free velocity field (To). Furthermore it is believed that 
the KPZ equation has the same large-scale behavior of 
the Kuramoto-Sivashinsky equation in 1 + 1 dimensions 
JlT| , while in higher dimensionality the situation is much 
less clear [Q. Nonetheless, in spite of the gigantic ef- 
fort devoted to the KPZ equation in the past decade, a 
complete understanding of its properties is still lacking. 

The aim of the present paper is to introduce a new in- 
verse approach to the KPZ equation. A previous attempt 
due to Lam and Sander |l3| was based on the standard 
Least Square (LS) reconstruction method. These authors 
used this approach directly on numerically simulated ex- 
perimental surfaces without a preventive test of the per- 
formance of the method itself. Lam and Shin |l4| subse- 
quently showed that the standard discretizations used in 
|13| was not adequate. We shall argue below that, even 
with the improvements given in [H], the classical iden- 
tification procedure devised in Ref. [ [i~3[ is not properly 
suited for Langevin dynamics since it is based on de- 
terministic equation ideas. By an explicit computation 
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using the LS technique applied to a 1 + 1 KPZ equa- 
tion, we shall review their method and point out what 
we consider its main deficiencies. 

We then go on to introduce a different approach based 
on the Fokker-Plank equation (FPE) associated to each 
Langevin equation jl5| , (J] . The advantage of this view- 
point is that one can construct deterministic relations 
among correlation functions which however still carry in- 
formations regarding the fluctuating nature of the orig- 
inal quantities. Those equations can then be easily an- 
alyzed within a least squares framework as in the LS 
method. 

The paper is organized as follows. In Sec. [ij], the 
KPZ equation is rapidly recalled along with its numeri- 
cal real space approximations in 1 + 1 dim ensions while 

Sec. " 



=28 id 6(t-t'). 



(4) 



the LS approach is reviewed in section III 



F 



tains the basic equations of our modified method which 
is then applied in Sec. [v]. Numerical results are then 
given in Sec. JVJ and some concluding remarks are pro- 
vided in Sec. [VI1| . More technical points are finally con- 
fined in the Appendices. Appendix shows why the least 
square method fails for sufficiently large noise amplitudes 
and Appendix |b] presents some results concerning renor- 
malized interfaces and their corresponding renormalized 
equations. 



II. INTERFACE DYNAMICS 

We consider a one-dimensional line of total length L 
and a surface of height h(x,t) at position x and time t. 
The continuum 1 + 1 KPZ equation then reads: 

d t h(x,t) = c + vd 2 x h{x,t) + ^[d x h(x,t)] 2 + V(x,t), (1) 

where rj{x,t) is an uncorrelated white noise 

(r)(x, t)rj{x', t')) = 2DS(x - x') 6{t - t'). (2) 

The average ( ) is taken on different realizations of the 
noise. In Eq.(0) and (||), c, v, A and D are coupling pa- 
rameters (c is often set to zero because of the invariance 
of Eq.(||) under rescaling h — ► h + ct). For A = 0, Eq.(|l]) 
reduces to the Edward- Wilinson (EW) equation Q which 
can be solved exactly. 

In writing Eq.(Q) cither a regularization in the corre- 
lation given in Eq. (|^) (such as for instance a spatially 
correlated noise ) or the introduction of a minimal length 
scale a is always tacitly assumed. In the latter case, one 
is then naturally led to consider a discretization of the 
continuum equation at a given cutoff length scale a. In 
that case, (a) the noise term rj{x,i) is discrctized 



with Si j the Kronecker symbol; and (b) Eq. (|l|) is written 
for a discrete variable hi(t) (i = 1, . . . , JV = L/a) with 
periodic boundary conditions 



SjL = c + u«Fr[h] + tg.F}[h] 



(•5) 



Here i/ c s = v/a 2 , A ff = A/a 2 , D c s = D/a. F"[h] and 
F*[h] are proper discretizations of the linear d 2 h and 
non-linear (d x h) 2 terms respectively. We note that the 
exact meaning of "proper discretization" has been the 
object of some investigations [^6|-^9|. 

In all practical applications, a further temporal dis- 
cretization [^o) is also performed on Eq. (|) 

hi(t + 6t) = hi(t) +5t(c+ Ves F?[h(t)} + F?[h(t)} 

where r*j is a Gaussian random generator of unit variance 
and 8t the discretization time step. 

In d = 1 + 1 it is known that the steady state solution 
P[h] for the probability distribution of the heights in the 
KPZ equation is identical to the EW stationary distri- 
bution due to a fluctuation-dissipation theorem ]l0[ . It 
was shown [Q, that the correct stationary discrete 
probability namely 



P[h] = Af^ 1 exp 



N 



(7) 



where Af 1 is a normalization factor, can be obtained by 
taking 



F?[h\ =h i+l + hi- x -2hi, 



(8) 



and 



(h i+ i - hj) 2 + Qn+i - hi)(hi - hi-i) + (h t - hi-if 



The standard choice F*[h] = 1/4 (h i+ i — /i;-i) 2 on the 
other hand fails to reproduce Eq. (ffl) and suffers of other 
problems as well |Q . A necessary (albeit not sufficient) 
condition for the identification with the continuum coun- 
terpart Eq. (Q), is clearly that the correct steady state 
(i.e. independent of A) is recovered. For this reason, we 
shall exploit for the new identification procedure as well 
as for the LS scheme, equations (|sj) and (^|) hereafter 
instead of the standard choice which was used in fll~3tl - 



III. THE LEAST SQUARES ERROR MODEL 
METHOD 



where 9i(t) is a random noise 



(3) Before introducing our method we first review the LS 

error method used in Ref. Jl3] |. We consider experimen- 
tal surfaces coarse-grained at length scale a described by 
the interface heights h° bs (t), (i = 1,...,N) which are 
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sampled M times i.e. at discrete times t = t k = kAt 
(k = 1, . . . , M + 1). Note that the sampling time At is 
the time interval between two experimental observations 
and it is clearly different from the discretization time St of 
Eq.(^J). For surfaces obtained by numerical simulations 
At is typically a multiple of St. We note that in Ref. 
Jul , the authors used At equal to St which is a rather 
particular case. 

For the sake of simplicity, we assume here that mea- 
surements are free from observational noise. It must be 
emphasized that, in the presence of measurement noise, 
our method performs a priori better then the LS scheme 
since it is based on spatial-averaged values which are less 
affected by errors on local height measurements. 

Our purpose is to determine the coefficients c, v, A, D 
at the given length scale a in Eq.(||). 

Let us first neglect the dynamical noise in Eq.(||). 
We then obtain a standard identification problem of the 
coupling parameters governing a deterministic non- linear 
equation which can be cast in the compact form: 



where we have defined 



dh 



(10) 



where in the present case p = 3 and \x\ = c, fi2 — v e B> 
M3 = A e ff, (0) and (g) are used for F?[h] and Ff[h] 
whereas Fl[h\ = 1. 

Optimal parameters are then determined by minimiz- 
ing a cost function J such as the sum-square difference 



J = 



1 



M N 



NM 



(11) 



fc=l i=l 



which quantifies the distance between experimental ob- 
servations h° hs (t k ) and equivalent reconstructed quanti- 
ties hf IC (tk)- The latter quantities are computed from 
Eq. (|10D for given parameters and are thus generally im- 
plicit functions of the parameters. However, if the sam- 
pling time At is small enough, then F^[h] are nearly 
constant between two measurements and the amplitudes 
/i? r (tfc+i) can be related to the parameters /i Q by 

hY ed {t k+1 ) = h° bs (t k ) + AtJ2 » a F?[h oh *(t k )] (12) 



M N 



fc=i t=i 



NM 



M N 

E£[- 

fe=l i=l 



h?*(t k+1 )-h? s (t k ) 



At 



Ff, 



(15) 



(16) 



and where functions F° are clearly expressed at 
hf> s (t k ),...,h$ s (t k ). 

This classical least squares method is an easy and natu- 
ral approach and it works fairly well in the absence of any 
noise. In the presence of noise however, it has its main 
drawback in the fact that it approximates time deriva- 
tives by finite differences. If the dynamics is governed 
by a deterministic equation and measurements are per- 
formed with a negligible observational noise, this simply 
imposes the choice of a sampling time much smaller than 
the characteristic or relaxation time of the process. 

Lam and Sander |l3| ] assumed that if the sampling time 
At is small enough, the above method could be extended 
to a Langevin equation (i.e. with dynamical noise). The 
amplitude of the noise can then be inferred from Eq. ( |Tl| ) 
when J is taken at the minimum values of the parameters, 
that is 



(17) 



However it was already observed in dynamical systems 
that even with pure measurement noise the above method 
can cause large errors. This is expected to be the case 
for dynamical noise as well. Two main reasons for this 
could be advocated. First if At is too large, the linear 
approximation ( p^ ) which explicitly relates the observed 
quantities breaks down. Because of the dynamical noise 
term, this happens a priori for shorter times intervals 
in a Langevin equation compared with its deterministic 
counterpart. Second, even in the favorable case in which 
At is small , such a method is efficient only if large sizes 
and small noise amplitudes are used. This is explained in 
Appendix ^ where a simple zero-dimensional case is ex- 
plicitly worked out with the method of Lam and Sander. 



In this case, the cost function J{{n}) itself becomes 
explicit and quadratic in the parameters. Optimal pa- 
rameters can thus be evaluated through a simple matrix 
inversion. Indeed the extremal value of J is inferred by 



dj 
dfi a 



=0. 



(13) 



The solution for the optimal parameters {/i*} is then 
given by a matrix equation: 



v 

= E] A oP ' B > 3 ' 

0=1 



(14) 



IV. STOCHASTIC APPROACH FOR MODEL 
IDENTIFICATION. 

We now turn to our method which is based on the sim- 
ple observation that all the information present in the 
Langevin equation (^) are also contained in the corre- 
sponding Fokker-Planck equation (FPE) |lj| : 

N „ 2 N „ 

d t P[h,t} = J2 D eB^joP[h,t]-J2^r( F i[h] P[h,t]), (18) 



where 
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Ei[h] = C + V eS F?[h} + -Aeffi^[/l], 



and 



JV 



p[h,t\ = (Udihi-hiit)) 



(19) 



(20) 



where the solution hi (t) is associated to a particular noise 
configuration 0j(t). 

In equation ( |l8|) the second term on the r.h.s. charac- 
terizes the deterministic behavior of the system whereas 
the first term contains stochastics effects. We derive a 
first general equation involving the parameters c and A. 
Using Eq. (fL8|), the time derivative of the ensemble av- 
erage of hi(t) can be easily shown to be: 



d(hi(t)) 



dt 



= J Vh Fi[h] P[h,t], 



(21) 



where Vh = ]X=i dhi. ^ we denote by g^(t) = 
77 Si (^(*)) the mean height at time t averaged over 
the noise, its time derivative can be written after some 



simple algebra as 
d 5 W(i) 



c • 



dt 6 
where we have defined 

N 



9i 



N 



(22) 



(23) 



in the variables Shi = hi — hi + \. Note that there are only 
N— 1 independent Shi variables due to periodic boundary 
conditions and to the fact that $^ i=1 Shi = 0. 

The above result prompts a convenient change of vari- 
ables from hi,..,hjy to Shi, ••> Shjy-i, h = l/NY^-^ hj 
followed by an integration over h. Physically this is re- 
lated to the fact that our system is infinitely degenerate 
with respect to the average height. Note that the station- 
ary probability Eq. (Q) is now Gaussian and well defined 
in the new variables Shi, . . . , The corresponding 

probability P[Sh] is the solution of a modified FPE: 



N-l 



d 2 



d t P[Sh,t] = 2D cS J2 5^2 P ^M 



N-l 



E^( Gi[5h]p[SKt] )- 2DeS ^ 



N-l 



d 2 



and 



l r 



Sh 2 _i - Shf +1 - Shi{Sh i+ i - Shi-i) 



(27) 



We are now in the position to derive our second basic 
result. 

Integrating Eq. (|24|) over all variables but (say) Shj, 
one gets, for the single variable probability 



p(Shj) = / VSh 3 P[Sh,t], 



(28) 



where the shorthand notation VShj = Yl&j=i dShi was 
again exploited, the following equation: 

d t p(Shj,t) = 2D eS -^p(Shj,t) - —niSh^t), (29) 

3 ^ 

where the non-local term ix{Shj,t) is defined as 

n(Shj,t) = / VShj Gj[Sh] P[Sh,t\. (30) 



The last step is to introduce the Fourier transform of 
p(Shj,t) which can be reckoned as a generating function 
for all moments of the distribution. Specifically on defin- 
ing 



Pj(Q,t) 



+ 00 



dShj e iqSh 'p(Shj,t), 



(31) 



we find a simple equation for the average p(q, t) over all 
sites oipj{q,i) : 



d t p{q,t) = -2D cS q 2 p{q,t) - iq7?(q,t), 



(32) 



in which 7c(q, t) is the Fourier transform of it{Shj,t) av- 
eraged over all sites. One can then expand Eq. ( |3^ ) 
in powers of q and obtain an infinite hierarchy (closure 
problem) in the correlation functions. The first two non- 
trivial orders (0{q 2 ) and 0(q 3 )) are 



dg£\t) 



dt 



= 4i/, 



err 



(2)/ 



-4D, 



ell, 



(33) 



and 



T, = — iv eS 

at 



g&\t)+g$(t)-2g&\t) 



+ 2 XcS 



g$i(t)-g&(t 



J5thp| ; gftwfi] have defined the following higher order corre- 
2 dShid5hi^ tion f un ctions 



where we have defined 

Gi = F. L — Fi + i 
with 

G\ = 5h i+1 + Shi-i - 2Sh 



v^Gl + -\ c ffGi, 



(25) 



(26) 



1 N 

sSW = ^E(««i(i) Sh l+l (t) 5h l+m {t)) , 



(35) 



i=l 



1 N 

= n E ( 5h *(t) Shi + l ® 5h ^) ^ +n (t)) . (36) 
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It is worth mentioning that A does not explicitly appear 
in equation (|35|). As one can explicitly check, this is 
a feature associated to the particular discretization Eq. 
(Q) and it would have not been the case had we used 
the standard discretization for F*[h]. This is clearly in 
turn related to the fact that the steady state probability 
distribution Eq. ([?]) is independent of A. We also note 
that in the 1 + 1 case we are considering, the explicit 
steady state solution of Eq.(|32) is known, and depends 
only on a single parameter ^ As a consequence, the 
steady version of Eq. (^) cannot be used here to identify 
v and D. In the 2 + 1 case where such a peculiar feature 
is not present, the stationary solution depends on A as 
well and parameters identification can exploit the steady 
analogue equation. 



V. PARAMETERS IDENTIFICATION 

Our aim is to implement an identification procedure 
which could be exploited on real experiments. For this 
reason, we assume that the experimental surface is con- 
stituted by a finite number of sites N with lattice spacing 
a (corresponding to a size L = Na), and it is measured 
during a finite time T a b s every sampling time At. Again 
At is a priori different from the discretization time 5t 
when the data is produced numerically (note that in real 
experiments St is not even defined!). We shall test the 
robustness and efficiency of the scheme with respect to 
the size L and sampling time At. 

Identification methods are often based on minimizing 
a cost function defined through dynamical constraints. 
This is clearly the case of the least square method as 
explained in section |lll|. H ere we derive dynamical con- 
straints using Eqns. (|2^) and (^) which contain in- 
formations of the original Langevin equation including 
mean values and fluctuations around mean values. The 
present identification is thus based on deterministic equa- 
tions. This constitutes a crucial difference with respect 
to the previous reconstruction method p3[ directly based 
on stochastic equations. Another important feature is 
that the observed quantities we use in our reconstruction 
scheme are dealing with averaged site values. Hence the 
fluctuations of all these terms, which derive from stochas- 
tic quantities, are reduced typically by a factor l/\/~N, 
and self-averaging is expected to be more effective. 

Let us derive the constraints we use. First, the to- 
tal observation time T b s is divided into q equal slices 



[Ti,T 2 ] 
tegrate 



l] with AT = Tj + \ — Tj. Let us in- 



2|) and ( |33| ) on each slice [Tj, Tj+i] : 



A3 



(i) 



If functions and integrals in Eqns. ( J37| ) and ( p8| ) are 
computed using experimental data, these discrete equa- 
tions provide 2q relations between the parameters to 
identify. From these constraints, two cost functions are 
built in a way already described in Sec.^ with p = 2. 
The corresponding 2x2 equations then yield c and A 
from one cost function and v and D from the other. 

We now explain how functions and integrals in Eqns. 
( |37| ) and (|3^) are obtained experimentally. Starting with 
a same initial surface e.g. a flat surface, we will grow 
the surface 1Z times. Because of the stochastic nature of 
the phenomenon, this produces 1Z different observations 
or realizations of the same process . Such a procedure, 
which can be performed very easily in real experiments, 
allows the computation, at sampling times t = tk = kAt 



Tj+i — Tj 



6 Tj+i - T , 



dt 



2gi 2 \t)+g[ 2 \t) 



(k = 1,...,M+ 1), of [ffWjexp, [^ 2) ] CXP and [g\ z Jcxp . 
Indeed, these quantities are the averages over 1Z differ- 
ent realizations of the spatial average height and corre- 
lations of the first neighbors. The number 7Z of realiza- 
tions need not be large: if the total number N of sites 
is sufficiently large, the experimental values are rather 
close to the corresponding theoretical predictions g^'(t), 

(21 (2) 

(7g (t) and g\ (t) . From these functions sampled every 
t = tk = kAt, integrals in Eqns. ( |37| ) and ( |38| ) can be 
efficiently evaluated for small sampling time At. In this 

case, the smooth functions [p^Jexpi [^o ]exp [g? ]exp can 
be approximated on the whole time interval [Ti,T 9+ i] 
by a standard curve fitting algorithm which gives as a 
by product the time integrals. This method does impose 
a constraint on the sampling time At. However, this 
constraint is substantially weaker with respect to that 
imposed by the LS method as we will show below. This 
is a considerable advantage of our new procedure. 

Two remarks are here in order. Firstly one expects the 
result to be independent on the number of slices q pro- 
vided that q satisfies the following two constraints. On 
the one hand q should be greater than 2 (since two param- 
eters are identified per cost functions) and on the other 
hand it should be less than M — -^j- so that AT cannot 
be less than the sampling time At. Secondly the identi- 
fication of v and D could be achieved by using Eq.(|32]) 
rather than Eq. (33). We shall see that in our approach 
the two equations yield virtually identical results. 



,(2)l 



(2) 



Tj+i — Ti 



T 3 + l - T j 



T 3 + l 



dt 



VI. RESULTS 

In order to test the potentiality of the different identi- 
fication methods, we produce experimental data by sim- 
ulating Eq. (^) with a standard Euler time integration al- 
gorithm with time step St = 0.01, lattice spacing a = 1, 
and3parameters v = D = 1 and A = 3. These are the 
same values used in Ref . [[l4| . The time step is expected 
to be sufficiently small for not causing instability prob- 
lems and the non-linear term A is big enough to be well 
^inside+he KPZ phase. We find interesting to repeat each 
lion few times (typically 5) to give an estimate of 
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the error bars to be associated to each parameter value 
(this was missing in previous works). 



A. LS Method 

Let us compute the parameters using the original 
LS method with the spatial and temporal discretization 
Eqns. (B) and (^). We exploit the same trick used in 
Ref. |jl4f] in which a KPZ surface of size 2L is obtained 
by a magnification of a fully relaxed surface of size L 
where the height are rescaled by a factor 2 Q , (a = 0.5) 
and linearly interpolated. The obtained surface is then 
relaxed to stationarity before a successive magnification 
is attempted. However, unlike Ref. Q where a single 
surface of size L = 32768 was computed, we consider 
L = 512,1024,2048,4096 and linearly extrapolate the 
results to the limit L — > oo. The calculation is repeated 
for increasing values of s = At /St in order to display 
the crucial weakness of the method as explained before. 
Fig. [j] depicts the results for the parameter v at finite L. 
Similar trends are present for A and D. The extrapolated 
values at L — ► oo are reported in Table [| The gradual 
decrease in the precision of the reconstructed parameters 
is apparent and it shows the loss of accuracy of the LS 
method as At increases as previously advertised. 

We also considered the LS when the reconstructed 
quantities are computed at time intervals AT which are 
multiple of the sampling time At. In fact this test was 
also carried out by the authors of Ref. |l3] (in their no- 
tations r = AT and At = St) and it will constitute a fur- 
ther source of comparison with our alternative stochastic 
method (see below). Even in this case there is a de- 
crease in the performance of the procedure as the ratio 
r = AT/ At increases, consistently with the results of 
Ref. jL4j. The corresponding extrapolated values are re- 
ported in Table y. 



B. Stochastic Approach 

For a more convenient comparison with the LS method, 
we use the same sizes and statistics (five different config- 
urations for each size). Our calculations are carried out 
in the transient rather than in the steady state and are 
therefore much less time consuming. Again the results 
are obtained for L = 512,1024,2048,4096 and linearly 
extrapolated to L — > oo. For a comparison with the pre- 
vious calculation, the outcomes of the parameter v at 
different size L are plotted in FigJ| for increasing values 
of the ratio s = At/ St, and the corresponding extrapo- 
lated values are reported in Table |l|. One can see that 
the parameter values are rather insensitive to the chang- 
ing the ratio s — At /St as expected. Next we check the 
performance of our method with respect to the increas- 
ing of the ratio r = AT/ At. This is reported in Table 



[V. As expected, our method outperforms the LS one in 
all situations. 

Since the LS method could in principle be carried out 
in transient rather than in steady state conditions, one 
might wonder how it would perform in this case. 

To this aim we recompute the parameters using the 
LS method under these conditions and find that the pre- 
dicted values are far off with respect to the exact ones. 
For instance for L = 4096 a typical run yields v ~ 0.36, 
A ~ 0.68 and D ~ 0.005 to be compared with a typical 
result obtained with our method v ~ 0.99, A ~ 2.98 and 
D ~ 1.01. 

As a final cross-check of our method, we recompute the 
parameters in the same situation as before but using Eq. 
( [§2l ) rather than Eq. (^) to extract v and D and find 
nearly identical values. 



C. Coarse-graining and KPZ real discretization. 

The application of this machinery to experimental sur- 
faces assumes that the system is described by a KPZ-likc 
dynamics. In this case, besides being able to address the 
issue of whether or not they belong to the KPZ univer- 
sality class, one would be able to provide a numerical 
estimates of the coupling parameters which are usually 
overlooked in studies focussing only on the universality 
class. 

Following Lam and Sander [ Of , we produce an inter- 
face based on the KPZ discretized model Eq. which 
is then smoothed by introducing the (discrete) Fourier 
transform of the heights 



h qn (t) 



N 
i=l 



(39) 



A coarse-graining surface at level a s = ba can then be 
achieved by simplying setting to zero all wavelength com- 
ponents hq n (t) with q > N s = L/a s and transforming 
back to real space. The obtained smoothed surface is 
then assumed to be governed by a KPZ equation (renor- 
malizability property). This new KPZ dynamics can then 
reconstructed along lines similar to those described above 
before a further time step is carried out on the original 
surface. 

With 6 = 2, AT/ At = 1 and sizes up to L = 8192 
averaged over 5 configurations as before, we find v = 
1.09 ±0.04, A = 3.27±0.05, and D = 0.88±0.03. Higher 
values of b result in poorer and poorer agreement with the 
expected values even with higher lattice sizes. The same 
feature is also present in the original LS procedure as we 
explicitly checked. In fact this is a general deficiency of 
the real space discretization as explained in Appendix 
the finite size difference has lost some renormalizability 
property of the original KPZ continuum equation. 
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VII. CONCLUSIONS 

In this paper, we discuss a method for extracting the 
coupling parameters from a non-linear Langevin equation 
starting from experimental surfaces representing succes- 
sive snapshots of the system. We apply this scheme to 
the KPZ equation in 1 + 1 dimensions (although it could 
be extended to any dimensions) and compare it with the 
previous approach of Ref . 113] , finding the following dif- 
ferences. First of all it does not require large sizes and 
it is well suited for a transient state. This is expected to 
be a considerable advantage, notably in numerical work, 
since the typical time required to reach a steady state 
increases as L z where z is the dynamical exponent (3/2 
in the 1 + 1 KPZ case). We have explicitly shown how 
the LS method which works rather well in the aforemen- 
tioned conditions, fails to provide sensible answers oth- 
erwise. Most importantly however is the fact that our 
approach is stable under the changing of the sampling 
time, unlike the LS method which is not. We stress the 
importance of this feature since in typical experimental 
situations, the sampling time is an externally tuned pa- 
rameters which has nothing to do with the evolution time 
of the system. We have discussed the reasons why this 
is so and provide an intuitive heuristic argument show- 
ing why the LS scheme is not expected to work under 
these more realistic conditions. Finally we implemented 
a coarse-graining procedure in order to be able to ap- 
ply our method to experimentally generated profiles. We 
showed that the agreement with the expected values is 
much poorer in the present case and we further argued 
that any real space based approach is doomed to run 
into this problem, the reason being that they have not 
a correct renormalization behavior under coarse-graining 
as explained in Appendix [§]. We have recently devised 
an alternative approach based on a Fourier-based scheme 
which avoids discretization problems. The results of this 
will be the subject of a forthcoming publication. 
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APPENDIX A: A SIMPLE SOLVABLE EXAMPLE 

This appendix shows, on a simple and solvable example 
which is a zero-dimensional analogue of Eq. ([!]) , that the 
method of least squares can be hampered by the presence 
of dynamical noise. Let us assume that the scalar variable 
X(t) is governed by the following Langevin equation: 



^=B(X) + »G(X) + r 1 (t), 



(Al) 



where B and G are prescribed functions, rj(t) is an un- 
correlated white noise 



( V (tW)) = 2D5(t-f). 



(A2) 



For simplicity, we assume that: (a) measurement noise 
is negligible, (b) the observed time serie X ohs (t k ) with 
t k = kAt (k = 1, . . . , M + 1) has been produced by the 



dynamical system (Al) with the value [i = 0. We ask 
whether the least square method is capable of identifying 
the correct coupling parameter fi = 0. 

From the one hand, the least square method first as- 
sumes that the da ta is produced by the deterministic 
counterpart of Eq.(Al) with an unkown parameter fi. In 
discrete times, this yields a "predicted" value given by 



X prca (t k+1 ) = X° os (t k ) + At[B(X oos (t k )) + ^G(X ot > s (t k ))\. (A3) 

On the other hand, the "observed" value is given, if the 
sampling time is s mall enough, by the discrete time coun- 
terpart of Eq.([Al]) with = : 

X ohs (t k+1 ) = X ohs (t k ) + At B(X obs (t k )) + r(t k )V2DAt, (A4) 

where r(t k ) is a gaussian random generator of unit vari- 
ance. 

Using both experimental observations X ohs (t k ) and re- 
constructed quantities X prcd (t k ), a cost function can be 
constructed: 



1 M 



(A5) 



k=l 



Using Eqns. (A3) and (A4), the cost function is readly 
rewritten as 



, m 

J =mJ2 [AtfiG(X ohs (t k )) - r(t k )V2DAi 



(A6) 



k=l 



The minimum value of J then satisfies the extremality 
condition 



which provides the value 

* = /2J gEfeli^fa)G(^° bs fa)) 
" V At ±J2iUG 2 (X° hs (tk)) 



Let us estimate the two sums appearing in Eq. 
Because of self-averaging, the denominator can be clearly 
rewritten as 



(A7) 



(A8) 



1 M 

E G 2 {X° hs {t k )) ps (G 2 ) , 



(A9) 



k=l 



the average being over the noise r\. The second sum in 
can be separated into two contributions 
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M 



M 



(**)] 



fe=l 



M 



fe=i 



$>(t fc ) G(A ohs (t fe )-(G) 



where using the Central Limit Theorem (CLT) 
have that 



we 



A/ 



(G) 



fe=i 



(All) 



In 



where ri is a random variable of unit variance 
Eq.( [A10| ), each term of the first sum of the r.h.s. is a 
random variable which is a product of two independent 
random variables of zero mean and variance respectively 
equal to 1 and a = (G 2 ) — (G) 2 . A similar argument 
based on the CLT applies to the evaluation of the first 
term of the r.h.s. since it is again a sum of M random 
variables with zero average which implies that: 



i£ Kfc) G(X° hs (t k )-(G) 



A I 



M 



k=l 



(A12) 



where r-i is a random variable of unit variance. Instead 
of the true value one then finds 



'2D 



Max[y/a, (G) 



VMM 



(G 2 



(A13) 



With non- vanishing amplitudes D, the noise hence causes 
errors for small At unless large statistics are considered. 



APPENDIX B: NON-RENORMALIZ ABILITY OF 
REAL SPACE DISCRETIZATIONS 

Let us consider a one-dimensional surface defined on a 
lattice of length L = N a (a being the lattice spacing and 
N being the total number of sites) which is identified by 
hi, ... , hpf at positions x\ = a, . . . , Xn — Na and having 
periodic boundary conditions. If we assume a KPZ dy- 
namics, then the equation of motion is given by Eq.(|5|) 
and its stationary state by Eq. ([?]) . Let us introduce the 
(discrete) Fourier transform h q so that 



N/2 



M*) = y E eiqnXi M*). 



(Bl) 



and conversely 



-N/2 



N 



By using the relation 



JV 

E< 

i=l 



i^qn-qm^xi 



N5 n 



(B3) 



it ig&a^$4c> show that the correct corresponding of Eq. (Q) 
jforl^c/starifeiai44llft)ribution is 



k=l 



P a [h] =A/'- 1 exp 



1 v 
2~DL 



N/2 
= -N/2 



(B4) 



In Eq.(B4) the continuum limit a — > is simply achieved 
by letting N — > oo. 

We now recall that the proper variables to be used 
in this context are the Shi = hi — hi+i whose Fourier 
transform Sh qn are related to the Fourier transform h qn 
of the heights by 



Sh gn =h qt Jl-e iq ' • 



(B5) 



We can exploit the periodic boundary conditions to ex- 
tend the TV — 1 sites to N (bearing however in mind that 
only N — 1 are independent) and use the the fact that 



N 



N/2 



i=l n=-N/2 



(B6) 



to rewrite the probability (B4) in the alternative form 



P a [h] =AA" 1 exp 



N/2 



2 DLa 2 ^ 

n=-N/2 



\K (i 



(B7) 



In the continuum limit a — > one can explicitly check 



that Eq.(B7) is equivalent to Eq.(B4) by expanding the 
term e lqna in powers of a and keeping th e lo west non- 
vanishing order 0(a). Nevertheless Eq.(B7) leads to 
problems when one tries to coarse grain the surface. In- 
deed suppose we now perform a smoothing of the lattice 
of a factor b to obtain a new lattice constant a s = ba and 
a decreased number of sites N s — N/b. This amounts 
to set to zero all modes from N s to N and thus the cor- 
rect corresponding stationary distribution is according to 



Pa s [h] =AC 1 exp 



1 



N,/2 



2 DLa 2 



E ( ] 



-NJ2 



(B8) 



On the other hand, had we started from the 
"smoothed" lattice and construct the differences 8h\ = 
h\ — h s i+b of the coarse-grained heights, we would have 
found 



( B2 ) P a , [h] = A/"" 1 exp 



1 



v 



N B /2 

2 DLa 2 ^ 

s n=-N„/2 



K„(i-e iq - ba )| 



(B9) 



which differs from the previous one. 

It is clear that this is a general problem of all dis- 
cretizations in real space. In the limit a — * in fact, one 



recovers from (B9) the correct expression (B4) 
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FIG. 1. The coupling parameter v for increasing lattice 
sizes L = 512, 1024, 2048, 4096, in the original steady-state 
LS method. All quantities are in dimensionless form. Error 
bars are of order of the symbol sizes and are consequently not 
displayed. Different curves refer to increasing values of the 
ratio s — At /St. The cross (X) indicates the exact value of 
the parameter v = 1. 
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FIG. 2. The coupling parameter v for increasing lattice 
sizes L = 512, 1024, 2048, 4096, as obtained from our re- 
construction method. Error bars are of order of the symbol 
sizes and are consequently not displayed. Different curves re- 
fer to increasing values of the ratio s = At /St. The cross (X) 
indicates the exact value of the parameter v = 1. 

TABLE I. Extrapolated values of coupling parameters 
^(oo), A(oo) and D(oo) as a function of s — At/St (see text) 
as computed from the original LS method (at steady state). 



Af /Av- 

L\lf 01 






UyOO j 


i 

5 

10 
20 


0.999 ± 0.003 
0.194 ± 0.004 
0.103 ±0.004 
0.049 ± 0.002 


2.980 ± 0.005 
0.595 ±0.015 
0.319 ±0.008 
0.139 ± 0.004 


1.000 ±0.001 
0.200 ±0.001 
0.100 ±0.001 
0.050 ± 0.001 


TABLE II. Extrapolated values of coupling parameters 
^(oo), A(oo) and D(oo) as a function of r = AT / At (see text) 
as computed from the original LS method(at steady state). 


AT/ At 


v(oo) 


A(oo) 


D{oo) 


1 

5 

10 
20 
50 


0.999 ± 0.003 
0.996 ± 0.002 
0.917 ±0.002 
0.855 ± 0.002 
0.650 ± 0.001 


2.980 ± 0.005 
2.817 ±0.004 
2.627 ± 0.006 
2.308 ± 0.006 
1.558 ± 0.006 


1.000 ±0.001 
0.928 ± 0.001 
0.853 ± 0.001 
0.757 ± 0.001 
0.653 ± 0.002 


TABLE III. Extrapolated values of coupling parameters 
v(oq), A(oo) and D(oo) as a function of s = At/St as com- 
puted from our stochastic approach in the transient state. 


At/St 


v(oo) 


A(oo) 


D{oo) 


1 

5 

10 
20 


1.009 ± 0.002 
1.008 ± 0.007 
1.035 ±0.011 
0.997 ±0.020 


3.047 ±0.016 
3.015 ± 0.006 
2.993 ±0.010 
3.001 ± 0.005 


1.026 ±0.001 
1.003 ± 0.007 
1.018 ±0.011 
0.978 ±0.010 


TABLE IV. Extrapolated values of coupling parameters 
v(oo), A(oo) and D(oo) as a function of r = AT / At as com- 
puted from our stochastic approach in the transient state. 


AT/ At 


v(oo) 


A(oo) 


D{oo) 


1 

5 

10 

20 

50 


1.009 ± 0.002 
1.003 ± 0.003 
1.003 ±0.001 
0.998 ± 0.002 
0.991 ± 0.003 


3.047 ±0.016 
3.005 ±0.010 
3.030 ± 0.007 
3.014 ± 0.009 
3.017 ±0.006 


1.026 ±0.001 
1.057 ± 0.003 
1.023 ±0.001 
1.016 ±0.001 
1.010 ±0.003 
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